Extended Gutzwiller Approximation for Inhomogeneous System 
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The generalization of the Gutzwiller approximation to inhomogeneous systems is considered, with 
extra spin-and-site-dependent fugacity factors included. It is found that the inclusion of fugacity 
factors reconciles the seemingly contradictory choices of Gutzwiller factors used in the literature. 
Moreover, from the derivation of the Gutzwiller factors, it is shown that the Gutzwiller approxima- 
tion breaks the rotational symmetry of the trial wavefunctions, and that different components of the 
spin-spin interaction need to be renormalized differently under the approximation. Various schemes 
to restore the rotational symmetry are discussed and are compared with results from variational 
Monte-Carlo calculations for the two-dimensional square-lattice antiferromagnet. Results along dif- 
ferent paths within the full parameter space, which corresponds to different choices of fugacity 
factors in the literature, are also compared. 



I. INTRODUCTION 

The t-J model in two dimensions has long been a fo- 
cus in condensed matter physics, for despite its seemingly 
simple appearance, it is believed to describe a variety of 
important strongly-correlated electronic systems, includ- 
ing the quantum antiferromagnets, the spin liquids, and 
the high-temperature superconductors fi The t-J model 
is described by the Hamiltonian: 



n 



(1) 



where P — Y\ji^^^3^^j^) the Gutzwiller projection op- 
erator, which accounts for the physical on-site Coulomb 
repulsion by explicitly prohibiting double occupation on 
any site. 

A common way to study the t-J model is to employ a 
variational approach, by introducing the Gutzwiller trial 
wavefunction: 



(2) 



where jV'o) is in general a Hartree-Fock type wavefunction 
that need not respect the double occupation constraint. 

With the Hamiltonian Ti. and the trial wavefunction 
lip), various expectations of the system can in principle 
be calculated numerically by the variational Monte Carlo 
(VMC) method.^ However, the VMC method IS compu- 
tationally costly and inefficient when long-range order of 
the system is sought, where the number of parameters in 
the trial wavefunction \ip) becomes large. 

A more practical way to perform these calculations is 
to employ the Gutzwiller approximation^ first introduced 
by Gutzwiller and subsequently clarified and extended 
by Ogawa et. al.^ Vollhardt)^ and Zhang et. alJ- In the 
approximation, an expectation with respect to {tp) is ap- 
proximated by multiplying the expectation with respect 
to \^po) by a factor that accounts for double occupation 
exclusion, i.e.. 



(Sj • Sj) w j)(S, • Sj)o 



where (Q) = (^|Q|?A)/(^|^) and (Q)o = 
('0o|'3|V'o)/('0o|'*/'o) for any operator Q. The expec- 
tation (Q) of the projected wavefunction lip) is said 
to be renormalized from the expectation {Q)o of the 
pre-projected wavefunction \ipo)^ 

Generally, the Gutzwiller factors gt{hj) and gj{i,j) 
are obtained by ignoring any non-combinatorial config- 
uration dependences of the expectation values. In the 
literature, for a homogeneous paramagnetic system, i.e., 
when (njt)o = (niT)o = (^ii)o = (?^^i)o for all sites i, j, 
the values of gt and gj are known to bei^i^ 



gt = 2S/{1 + S) 
5j = 4/(1 + 5)2 



(4) 



where 6 = 1— (ftij.)o ^ ("-11)0- Alternatively, Eq. |3]can 
also be derived from a functional integral approachjiSi 

However, the generalization of the two factors in an 
inhomogeneous case is not obvious. Some author a^^d^d^ 
simply take the results from Eq. [?] and interpret S as site 
dependent. Hence, they obtain: 



gS,j) = ^AS,s,/{i + s,)ii + Sj) 

gjit,j)=A/il + S,){l + S^) 



(5) 



where 6i = 1 — {ni)o is the local hole density at site-i. 

In contrast, for an anti- ferromagnetic trial wavefunc- 
tion in a square lattice. Can et. ali^ obtained: 



gtiA,B) = n{l - n)/{n- 2n+n_) 
gj{A,B)=n^/{n-2n+n_f 



(6) 



where A, B are labels for sublattices, and which = 
("AT>o = ("-si)o, = {fiAilo = ("-ST>o, n = n+ + n^. 
It should be noted that the gt they obtained is identical 
to that of Ogawa et. al.^ who did not derive gj. 

Another suggestion for gt and gj is given by Wang 
et. al.f^ who derived their results from grand-canonical 
wavefunctions and considered a generalized Gutzwiller 
(3) projection operator P' = WjV^^il — fij^fiji), where yj 
are local fugacities determined by the condition that 
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{nj)o = {fij). For the t-J model, the factors they ob- 
tained are of the form gt{i, j, a) = -Jgt{i, <^)y/gt{j, c) and 
9j{hj) = y/9j{i)^/9j{j), where: 

N I [l - n,g){l - n.,^ - n^i){ni^ + nil) 

Jgjyi) ~ ■ 

with ni(j — {nia-)o- These equations reduce to Eq. |3] in 
the homogeneous paramagnetic case, and to Eq. [S] in 
the square-lattice antifcrromagnet. However, in a ho- 
mogeneous partial ferromagnet, Eq. [7]is in disagreement 
with the results from Zhang et. alX, which are given by 
9t{cr) = (l-n|-nj/(l-n^) and = l/(l-7i|)(l-7ix). 
It should be noted that this gt obtained by Zhang et. al. is 
identical to that of Ogawa et. al.^ who did not derive gj. 

The brief survey above indicates that there has been 
confusion in what the appropriate Gutzwillcr factors 
should be for inhomogencous systems. The purpose of 
this paper is to clarify such confusion by deriving the 
Gutzwiller factors by appropriately generalizing Ogawa's 
original approach^ and to investigate the accuracy of 
the resulting Gutzwiller approximation in inhomoge- 
ncous systems by comparing with VMC results in two- 
dimensional square-lattice antiferromagnct (SLAF). 

Our derivation, which is presented in the appendix, 
is based on configuration counting on canonical (i.e., 
particle-number eigenstate) wavefunctions. This is a 
more natural choice for deriving the Gutzwiller factors, 
since the Gutzwiller approximation amounts to neglect- 
ing quantum correlations in configurations, but retaining 
combinatorial ones. This combinatorial dependence is 
more clear in canonical wavefunctions as compared to 
grand-canonical ones (readers may want to compare our 
derivation with that of Ref. ITsl ). 

The key insight that resolves the confusion is that 
the single-particle density will in general be modified by 
the Gutzwiller projection, and that the relation between 
the projected density (fii^) and the prc-projected density 
(nicr)o can be adjusted by introducing local fugacity fac- 
tors Uia- that in general depend on both site and spin. 
This spin-and-site dependent fugacity factor was first 
introduced by Gebhardi^ when deriving the Gutzwiller 
factors using a diagrammatic approach. Similar factors 
also appear in works related to Gutzwiller projected su- 
perconducting stateSf^iiiii^ where the fugacity factor is 
introduced to keep the mean number of particles un- 
changed. It should be remarked that the spin dependence 
of 2/icr is essential, since otherwise one can merely main- 
tain {hi) — (rii)o, as in Ref. [isl, but not (hi^) = {hia)o, 
and hence the local magnetization {rhi) = {hi^) — {hii) 
may be modified by the projection. When making the 
physical argument that the Gutzwiller factor is given 
by the probability for the process to occur in the pro- 
jected state divided by the corresponding probability in 
the prc-projected stated the renormalization of single- 
particle density must be taken into account. 



In the following sections, we shall find out that both 
Eq. [5] and Eq. [5] follow from the derivation, but corre- 
spond to two different implicit choices of local fugacity 
factors. Specifically, the former results from demanding 
that {njcy)o = (fija) for all sites j and for all spin a, while 
the latter results from setting yj^a- — 1 for all j and a. In 
the particular case of homogeneous paramagnet, our re- 
sults reduce to Eiq.\^regardless of the value of y, which by 
symmetry is spin and site independent. This is expected 
for a canonical wavefunction and is in contrast with the 
derivation by Wang et. al.^ which demands a particular 
value of y for the equations to work out. An understand- 
ing of this implicit choice of fugacities is particularly im- 
portant if we want to compare results from Gutzwiller 
approximation to that from VMC calculations, since we 
need to make sure that we are comparing expectations 
with respect to the same wavefunction. 

We shall also discover subtleties associated with the 
definition of gj in scction|TTl In particular, the Gutzwillcr 
factor gjz for the z-componcnt of the spin-spin inter- 
action is in general different from the corresponding 
Gutzwiller factor gjxy of the x- and y-component. The 
physical origin of this difference is that in a fixed spin 
basis, the Pauli exclusion principle posts a restriction on 
the legitimate configurations in the pre-projected wave- 
function for the exchange of opposite spins, while such 
restriction is absent for exchange of same type of spins. 
Consequently, the Gutzwillcr approximation breaks ro- 
tational symmetry even in the homogeneous case, and 
Eq. m can only be reached by enforcing rotational sym- 
metry beyond the probability ratio argument. 

It should be noted that the configuration counting 
approach we used here is in general different from the 
l/d-expansion developed by Metzner, VoUhardt and 
GebhardJ^ii^ It is however the configuration counting 
approach that corresponds to the physical intuition that 
the Gutzwiller factors arc obtained by dividing the prob- 
ability that a process would occur in the projected wave- 
function by the corresponding probability in the pre- 
projected wavefunction (c.f. Sect, im and Appendix). In 
the case of SLAF, the formula for the various gj given by 
the two formulations are different, even when the d oo 
limit is taken in the l/d-expansion. However, the nu- 
merics agree qualitatively, and the numerical differences 
between the two approaches is of the same order of mag- 
nitude as the error between each approximation and the 
VMC result (c.f. Fig.[6]in Scct.Hnj. 



II. INTUITIVE ARGUMENTS FOR 
GUTZWILLER FACTORS 

In the homogeneous paramagnetic case, an intuitive 
picture of the Gutzwiller approximation is provided by 
Zhang et. al.-^ whose central result is that the Gutzwiller 
factors are given by the probability for the physical pro- 
cess under consideration to occur in the projected wave- 
function divided by the probability for such process to oc- 
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Physical process 



Probability in ip 
(before process) 
Probability in -0 
(after process) 



TABLE I: Probability for various physical processes to occur in hopping 



Probability in ?/;o (ni^)o{hii)o {nii}o{nii)o {ni^)o{l - mijo (rii|)o(l — ?iii)o 

(before process) x{l - nj^)o{nji}o x (1 - rij|)o(l - nji)o x{l - hj^)o{hji}o x{l - hj-i)o{l - hji)o 

Probability in i/;o (1 - ?iiT)o("-a)o (1 - n,|)o(?iii)o {1 - nit}o{l - nii)o (1 - n,|)o(l - ni|)o 

(after process) x ("-iT)o(^^ji)o x (nj|)o(l - nj|)o x ('^jT)o(nj|)o x (7ij|)o(l — nj|)o 



{fii^){'^-njy -hji) 

(1 - Mil - nii){nj^) 



cur in the pre-projected wavefunction. A careful deriva- 
tion of the Gutzwiller factors from the approach of Ogawa 
et. al.^ which we relegate to the appendix, confirms this 
intuitive picture even in the inhomogeneous case, with 
the exception that one can no longer assume the pro- 
jected density {'hia-) to equal to the pre-projected density 
{nia)o- In this section, we shall review the argument by 
Zhang et. ai, from which we shall see that gj, ^ gjxy 
We shall also supplement the result of Zhang et. al. by 
providing an intuitive picture of how {hia) is computed 
in the inhomogeneous case. For notational simplicity, 
henceforth we denote {hi^) as pia- and {hia-)o as rii^. 

We shall work with the generalized Gutzwiller pro- 
jector P' = Uj Vj^^ yji^ - ^iT"-ji)' where y^^ are 
positive parameters that depend on both site and 
spin. To avoid the complications from particle number 
renormalization^Siiiii^ throughout this paper lipo) is as- 
sumed to be a spin-definite canonical wavefunction (i.e., 
fixed number of up-spin and down-spin particles). 

First consider the hopping {cj^Cia). For concreteness 
take a ='\. In the pre-projected wavefunction, before 
hopping an up-spin must reside on site-i, and Pauli exclu- 
sion principle demands no up-spin on site-j. Since there 
is no occupation constraints, one or zero down-spin can 
reside on each site. After hopping, the down-spins are un- 
affected while the up-spin moves to site-j, leaving site-i 
with no up-spins. There are thus four legitimate config- 
urations in the pre-projected wavefunction. However, as 
a result of occupation constraint, only one of these four 
configurations is allowed in the projected wavefunction, 
namely the one which both sites contain no down-spins 
(see Table. U for illustration). Summing over the proba- 
bilities of legitimate configurations, taking the ratio, and 
include an overall square root, we have: 



r(l - rija) 



(8) 



This result is also produced and discussed in Ref . d. 

Next consider the spin-spin interaction. In a fixed spin 
basis, the spin-spin interaction consists of four distinct 
types of physical processes, corresponding to the different 
ways of expanding Si ■ Sj in terms of c and c^, and various 



(a) 



(b) 



(c) 



FIG. 1; Legitimate configurations in spin-spin interaction for 
the pre-projected wavefunction for various processes. Pro- 
cesses outside the rectangle in broken lines are disallowed in 
the projected wavefunction. (a) The process (c]|Cj|cJj^Ci|) 
(type 1). (b) The process (c]|Cj|cJ|Ci|) (type 2). (c) The 
processes in (cJ|Cj|cJ|Ci|) (types 3 and type 4). 



ways of Wick-contracting the four-fermion expectations. 
These processes are: 

1. Exchange of an up-spin with a down-spin (^ 

2. Counting of one up-spin and one down-spin ('^ 

3. Exchange of two spins of the same species 

{Cj^Cj-^C-^Ci^)) 

4. Counting of two spins of the same species (~ 



Note that the first type of processes contributes to the x- 
and y- component of the spin-spin interaction, while the 
remaining three contribute to the z- component. 

For processes of type 1, Pauli exclusion principle posts 
a strong restriction on the legitimate configurations in 
the pre-projected wavefunction. In the specific example of 
(cl-^CiiCj^Cji) , it demands site-i to have no down-spin and 
site-j to nave no up-spin. Consequently, the pre-projected 
wavefunction has only one legitimate configuration for 
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this process. In comparison, processes of types 2-4 all 
have four legitimate configurations in the pre-projected 
wavefunction. In other words, the occupation constraint 
is automatically demanded by processes of type 1, while 
has to be imposed additionally for process of types 2-4. 
See Fig. [T] for illustration. 

Consequently, the probability for processes of type 1 
to occur in the prc-projcctcd wavefunction has an extra 
factor of (1 — njcr)(l — njg){l — nig){l — nj„) when com- 
pared with the other types of processes. Repeating the 
exercise of calculating the overall probability for a given 
type of process and group together various terms, the 
renormalization factors for the spin-spin interaction can 
be summarized as: 

(Si • Sj) = gjxy{Si+Sj- + Si-Sj+)n + gjz{Siz)o{Sjz)o 
, g,/T ,4.1,. 41 , gji ,4 . , 4 , 



gjl = gjxy "by hand." However, since in the homo- 
geneous case {Siz)o = 0, we cannot similarly argue for 
setting gjz = gjxy 

To compute the Gutzwiller factors, we still need to re- 
late pi„ to riiCT- Again we relegate the more mathematical 
derivation to the appendix and focus here on physical in- 
tuitions. Given any pre-projected spin-definite canonical 
wavefunction l^/'o), we may expand in the configura- 
tion basis. Those configurations with sites occupied by 
both up-spin and down-spin will be projected away by the 
Gutzwiller projector P' . Since the Gutzwiller approx- 
imation conceptually amounts to neglecting correlation 
between different sites, each configuration that survives 
after projection should be assigned a classical weight W , 
which is a product of weighing factor w{i) on each site i. 
The appropriate expression for w{i) turns out to be: 



(9) w{i) 



where, 



gjxy 



gjz 



gj] 



gji 



P^lPjl 



(10) 
(11) 
(12) 

(13) 



Note that the last three terms in Eq. [9] arc all contribu- 
tions from the z-component of the spin-spin interaction, 
while the first term is the contribution from the x- and 
y-component. 

The four Gutzwiller factors above are in general un- 
equal. In the specific case of a homogeneous param- 
agnet, pia- ^ ni„ and hence g,jz = g.j^ = 5Ji = 1 
while gjxy = 4/(1 + 5)^. However, since the Gutzwiller 
wavefunction for a homogeneous paramagnet is rotation- 
ally invariant, physical reasoning demands {SizSjz) ~ 
gj\{c]^c,^)o{cj^cl^)o/4: + gji{cl^c,i)o{cjicl^)o/4: to equal 

to (SixSjx) ~ gjxy{Si+Sj- + Si-Sj+)o/2. As the pre- 
projected wavefunction is also rotationally invariant, we 
must also have {SizSjz)o = {SixSjx)o- Since {SizSjz)o 
is simply (SizSjz) with all Gutzwiller factors set to 1, 
and similarly for {SixSjx)o, one should expect gj-^ ~ 
gjl = gjxy in the homogeneous case. In other words, 
the Gutzwiller approximation derived from configuration 
counting breaks rotational symmetry. In retrospect this 
is not all surprising — the counting of configurations and 
the computation of probabilities are all done assuming a 
specific spin basis, and hence there is no a priori reason 
to expect rotational symmetry to be preserved. 

When applying the Gutzwiller approximation to near- 
homogeneous system, it may be desirable to set gj^ ~ 



nia-{l — nia)yf„ for site occupied by spin cr 
(1 — Jiif )(1 — n,;|) for empty site 



(14) 

The above expression for w{i) make conceptual sense: 
nia(l — riig) and (1 — rii^){l — riii) are respectively the 
probabilities for finding a cr-spin and an empty site at 
site-i in the pre-projected wavefunction, and the factor 
of yf^ comes from the factor of y"^' in both the bra and 

theketof (V-|V'> = (Vo|i^"|V'o). 

After assigning each configuration with a classical 
weight, the projected density pi^ is simply given by the 
weighed average occupation of site-i by a cr-spin. 

To derive an explicit set of equations that relate pi^ 
to riia, assume that the system consists of k sublattices 
(labeled by /), such that on each sublattice riic ~ nj„ 
and yiij = yj„. Let N be the number of lattice site and 
M = N / K. Then, for a configuration with a/^ up-spins 
and ail down-spins on sublattice /, the classical weight 
is: 



W{{ai,}) = X{{{l-nji)(l-nii)) 



^{y],ni,{l-nn))^'\yj,nn{l-ni,)f'' 

(15) 

Moreover, simple combinatorics shows that the number 
of such configuration is given by: 



c(W4)=n 



M! 



(16) 



For a fixed site i in sublattice P, a fraction of ap^/M 
out of the C({a/o-}) configurations contain a cr-spin on 
site-i. Hence, the weighed average is given by: 



Pier = 



E(a,,}(ap./M)C({a,4)iy({aj4) 
T.{a,.)C{{ai,})W{{ai,}) 



(17) 



In the thermodynamic limit where N 00, the func- 
tion F{{ajcr}) = C{{aicr})W{{aicr}) is sharply peaked 
and hence the sum J2{ai„} replaced by the sin- 

gle term in which F({aja-}) attains maximum under the 
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constraints '^jai^ = (V'o|'T'jct|'0o) = No-. The maxi- 
mizing o/g. can be found by standard Lagrange multiplier 
technique. Upon simplification, it is easy to check that 
Pia- can be solved through the following set of equations 
for a given choice of yta : 



Pla 1 - nia 



(18) 



here Ao- are Lagrange multiplier to be determined from 
the set of equations. Note that multiplying each by 
the same constant will not affect the equations relating 
pic to n/o-, and similarly for multiplying the same con- 
stant to each y/x*^ 

With the choice yia = 1/(1 — nia-)/{l — nj-^ — n/^), we 
have {hiu) = {nia-)o within the Gutzwiller approximation. 
With this implicit choice of y/o- then, we get: 



(l-"»T~'^a)(l~"jT~"ji) 
(1 - njo-)(l - Uja) 

1 



gjz = gji = Ji = 1 



(19) 

(20) 
(21) 



If we assume n^i = = {hi)o/2, we recover the gt and 
gjxy quoted in Eq. [3 However, we now see that to be 
consistent with the approximation scheme, one needs to 
set gj^ = gj^ = gji = 1. 

Note that the geometry and dimensionality of the un- 
derlying lattice have never entered the above discussion. 
This is a consequence of the approximations made in the 
configuration counting approach, where all quantum cor- 
relations except combinatorial ones are neglected. Com- 
binatorial factors in general do not depend on the de- 
tailed geometry and dimensionality of the underlying lat- 
tice. 

Since the results we obtained is independent of dimen- 
sionality, one may wonder if it is equivalent to the d 00 
limit of the l/d-expansion developed by Metzner, VoU- 
hardt and Gebhardj ^^'^^ For the renormalization factor 
gt, the approximations made by the two approaches are 
essentially the same, and the resulting factors are numer- 
ically equali^ For the renormalization factor gj, which in- 
volves a four-fermion expectation, the agreement is less 
spectacular. Although the derivation of Gcbhardi^ also 
indicates that gjz, gjxy, and gja^ should in general be 
different, in the case of SLAF the formula for gj^ from 
Gebhard's derivation disagrees with both gj^ and gjxy 
we gave in Eg . [iTIl [T51 It should however be noted that, at 
least in the case of SLAF, the numerical energy estimates 
from the two formulations agree qualitatively, and the nu- 
merical differences between the two approaches is of the 
same order of magnitude as the error between each ap- 
proximation and the VMC result (c.f. Fig.[S]in Sect. lIIip . 



III. THE CASE OF TWO-DIMENSIONAL SLAF 

To determine the accuracy of the Gutzwiller approx- 
imation and various modification schemes of it that re- 
store rotational invariance for the paramagnetic state, 
various expectation values in the two-dimensional SLAF 
are computed in these approximation schemes and are 
compared with the results obtained from VMC. In the 
following, we shall first treat yia as an additional vari- 
ational parameter and consider the entire parameter 
space, from which we shall discover that there is a vast re- 
gion in which the Gutzwiller approximation is erroneous. 
However, we shall also see that the expectation values de- 
pend strongly on physical density {fii^) and only weakly 
on how such density is obtained from the parameters. 
Consequently, there are paths within the full parameter 
space in which the full range of (hicr) can be explored 
and which the erroneous region can be avoided. We shall 
then focus on some of these paths and explicitly consider 
the different modification schemes. 

For fixed hole density 6 = ~ '^U ~ "-it)/^' ^^'^ 

prc-projectcd SLAF state is uniquely characterized by 
the antiferromagnetic order parameter A. Explicitly, the 
pre-projected wavefunction \ip) is given by: 

l^>= n K5L+sign(a)«k4+Q,J|0) (22) 
where Q = (tt, tt), ek = — 2(cos + cos ky), 



(23) 



and which the Fermi energy e^? is determined by the 
dopping S via Y^^{nicr)o = 1-5. 

This prc-projccted wavefunction is invariant under the 
simultaneous exchange of sublattice and spin indices. 
To preserve this symmetry after projection, we demand 
VAcr = VBa, where A, B label sublattices. Since multiply- 
ing all fugacity factors by a constant amounts only to an 
overall normalization of the resulting projected wavefunc- 
tion, the projection P' = Y\- Vj^^ l/jl^ (1 ~ ^jT^-iJ.) ^^^^ 
case is uniquely characterized by = yB-^/yAi- Hence, 
the wavefunction \tp) to be considered in this section is 
uniquely characterized by the physically adjustable S and 
the two trial- wavefunction parameters (j/r, A). Further- 
more, since the mapping {yr,A) 1— > {y~^,—A) amounts 
to inverting the roles of sublattices A and B, we need to 
consider only the parameter space in which A > 0. 

Unless otherwise stated, the explicit data we present 
are at dopping S = 0.025, and the VMC calculations are 
performed in a lattice consisting of 8 x 10 sites, with 
periodic boundary conditions. We shall first make the 
comparisons over the whole parameter space (y^, A) and 
then focus on specific paths within the parameter space. 
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FIG. 2: (Color online) Staggered magnetization (m) in (a) 
the VMC calculation and (b) the Gutzwiller approximation. 
The broken curve in (b) indicates the parameter subspace in 
which pia = riia within the approximation. The VMC data 
is interpolated from a grid of interval 0.1 in y,. and 0.05 in 
A. The typical error in (m) in the VMC calculation is about 
0.01. 



For two-dimensional SLAF, Eg. [T51 gives: 

n^(nj^ -f n„)(l — n_) 



PA\ = PBi 



PB] = PAi 



(24) 



where n+ — uai — nsi and n_ = uai — tibi as in 
Eq. [HI It should be remarked that by setting yr = 1 and 
plugging into Eq. [TOlfTSl we find gjz = gjxy and recover 
the results in Eq. [6]for gt and g,j {— gjz = gjxy), while 
5JT = 9.JI = (1 - ~ n_)gj. 

From Eq.[24]the staggered magnetization (m) = pAi — 
PAi can immediately be evaluated. The results and 
the comparison with VMC results are shown in Fig. [21 
Observe that constant magnetization contours from the 
Gutzwiller approximation are generally flatter than that 
from the VMC calculation. 



Next we consider the hopping expectation 



for 



nearest- neighbor sites i, j, which also gives (7i) in the t- 
model where J = 0. From the properties of the wavcfunc- 
tion, this expectation is independent of a and the choice 
of particular sites. The results are shown in Fig. [31 It is 
important to note that while the hopping expectation cal- 
culated from VMC shows a shallow global minimum near 
(yr, A) = (1,0), the hopping expectation calculated from 
the Gutzwiller approximation keeps decreasing along the 
direction in which both y^ and A increase to maintain 
an approximately unchanged small staggered magnetiza- 
tion. This decrease is unphysical and indicates a sys- 
tematic error in the Gutzwiller approximation (as can 
be seen in Fig[3Jc)). The origin of this high-error re- 
gion in parameter space can be understood as follows: 
since the Gutzwiller approximation is based on neglect- 
ing non-combinatorial configuration dependences of ex- 
pectation values, the effect of non-homogeneity caused 
by the fugacity factor yi^r (which is purely combinato- 
rial) can be accurately accounted for in the approxima- 
tion, while the effect of non-homogeneity caused by the 
parameters within the pre-projected wavefunction lipo) 
such as A cannot. The high-error region corresponds to 
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FIG. 3: (Color online) Nearest neighbor hopping expectation 
(cJ^Cjct) in (a) the VMC calculation and (b) the Gutzwiller 
approximation, together with (c) the magnitude of difference 
between the two. The VMC data is interpolated from a grid of 
interval 0.1 in yr and 0.05 in A. The typical error in (cl^Cja) 
in the VMC calculation is about 0.002. 



trial wavefunctions in which each of yi^ and A alone cre- 
ates a large non-homogeneity, but which the two cancel 
out each other to produce an almost homogeneous state. 
In these states the effect of each of yi„ and A on {c\^Cja) 
is large, but the former is estimated accurately while the 
latter is estimated erroneously, thus producing a large 
overall error. Note also that the error in the Gutzwiller 
approximation becomes small when A 3> j/r- This is not 
surprising as the wavefunctions become almost classical 
antiferromagnets in such limit. 

Next we consider the expectation of the Hamiltonian 
7i as defined in Eq. [H where we take the conventional ra- 
tio t/ J = 3 and without the loss of generality set J = 1 . 
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-0.25 




FIG. 4: (Color online) Energy Expectation (7i), where Ti. 
is defined as in Eq. [T] with t=3 and J=l, in (a) the VMC 
calculation and (b) the Gutzwiller approximation, together 
with (c) the magnitude of difference between the two. The 
VMC data is interpolated from a grid of interval 0.1 in yr and 
0.05 in A. The typical error in {Ti) in the VMC calculation 
is about 0.004. The white cross in (a) indicates the location 
of global minimum in the VMC data set. 



The results are shown in Fig. 31 where the Gutzwiller 
approximation of the J-term is implemented by Eq. llOj - 
[T51 Again we see that whereas the expectation ealculated 
from VMC shows a shallow global minimum, the expecta- 
tion calculated from the Gutzwiller approximation keeps 
decreasing along a direction where both yr and A in- 
crease to maintain an approximately unchanged small 
staggered magnetization, which again indicates systemic 
error in the Gutzwiller approximation. By splitting the 
contribution between the t-term and the J-term, it can be 
verified that this systemic error exists in the Gutzwiller 



4aCja) and {H) 



approximation of both {c\^Cj^) and (S^ • Sj). Moreover, 
whether we replace gjf and gj^ by gjxy does not affect 
the existence of this systemic error. 

Observe that in the VMC results, 
along a constant staggered magnetization contour are ap- 
proximately constant. It is thus sensible to speak of the 
various expectations as functions of {fiia), and to com- 
pute these quantities along a particular path within the 
whole parameter space, as previous authors have implic- 
itly doneiiiii^ii^ii^ Although the Gutzwiller approxima- 
tion is inaccurate in certain regions of the whole parame- 
ter space, it may be acceptable along some specific paths. 
In particular, the inaccurate region can be avoided by 
either paths with a fixed yr close to 1 or by the path 
where pi^ = ni„ within the Gutzwiller approximation. 
We shall now compare results from the Gutzwiller ap- 
proximation and the VMC calculation along these paths, 
and explicitly consider various schemes that repair rota- 
tional symmetry. For paths of fixed yr , we consider the 
following three schemes: 

(z,z) gjz, gjxy, gji and as obtained from Eq. 1241 and 
Eq. MM 

(z,xy) gjz and gjxy as in scheme (z, z), but with gj| = 

gji = gjxy 

(xy,xy) gjxy as in scheme (z, z), but with pj| ~ gji = 

g.jz = gjxy 

For the path along which pia = within the Gutzwiller 
approximation, we have one additional scheme: 

diagrammatic Gutzwiller factors obtained from dia- 
grammatic approach;— where gjz and gjxy agree 
with scheme (z,z), and which gj^ and gji are given 
by: 



5T ^91 



(l-n++n_)(l 



n-) 



(l-n+)(l-n_) 



(25) 



Since we are interested in whether the Gutzwiller ap- 
proximation can arrive at sensible physical results, we 
shall compare between the Gutzwiller approximated ex- 
pectations as functions of (projected) staggered magneti- 
zation within the approximation and the VMC computed 
expectations as function of the VMC computed staggered 
magnetization. 

First we consider the nearest-neighbor hopping expec- 
tation (cl^Cja-), which relates trivially to the energy ex- 
pectation in the t-model (i.e. the t-J model with J set 
to 0). The plots of (cl^Cja-) as function of (m) along 
the path where yr = 1 and pia- = are shown in 
Fig. [S] The most striking observation from the plots 
comes from the Gutzwiller approximated result along 
Pia = n.ic, which exhibits a minimum at (to) « 0.70, 
contrary to the results indicated by the VMC calculation. 
This false minimum is likely to be the remnant of the ef- 
fect that produces the large error region, since the curve 
where pi^j = is in close proximity with the large error 
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FIG. 5: (Color online) {cl^cjcr) as function of (?h) from the 
Gutzwiller approximation (unbroken purple line) and from 
the VMC calculation (blue filled dots), for (a) path along 
which yr — 1 and (b) path along which pi^. — riia. 



region near (rh) w 0.70 (c.f. fig. [2] and [3]). Moreover, ob- 
serve that the variation of (cl^Cja) as a function of (m) for 
small values of (m) in the VMC results is much smaller 
than the error between the VMC and the Gutzwiller re- 
sults. Thus a small error in the Gutzwiller approxima- 
tion for 



can lead to an incorrect prediction of 
the staggered magnetization at the t-model minimum- 
energy state. In other words, the systematic error in the 
Gutzwiller approximation may misled one into believing 
that an antifcrromagnctic state is stabilized even in the t- 
modcl (sec also the remark at the end of this section). In 
light of this, the result that spontaneous spin and charge 
ordering occurs in the triangular-lattice t-model2S should 
be viewed with some caution. 

However, it should be noted that even without com- 
parison with the VMC results, this inaccuracy could have 
at least been suspected, since the shape of {c\^Cja) as a 
function of (m) and the corresponding location of min- 
ima disagree between results from Gutzwiller approxi- 
mation along the two different paths in the parameter 
space. In other words, checking the results for different 
paths within the parameter space may provide a consis- 
tency check within the Gutzwiller approximation. 

Next consider the energy expectation {TL). The plots 
of {Ti) as function of {rh) along the paths where yr = 1, 
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FIG. 6: (Color online) (H) as function of (rh) from the 
Gutzwiller approximation (with various schemes described in 
text) and from the VMC calculation, for (a) path along which 
j/r = 1, (b) path along which yr — 1.1, and (c) path along 
which piCT = Uicr. The inset of (c) plots the curve for scheme 
(xy,xy) in its full range. 



yr — 1.1 and pi^ ~ are shown in Fig. [6l In this 
case, the VMC results consistently indicate a magnetized 
lowest-energy state around (m) « 0.78. Moreover, the 
shapes of curves from the Gutzwiller approximation gen- 
erally resemble the shapes of the VMC curves, with the 
notable exception of the (xy,xy) scheme along pi„ = . 
Recall that the rotational invariance for the paramag- 
netic state suggests merely that gj^ = gji = gjxy, but 
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bears no conclusion for gjz, since {Siz) vanishes for such 
state. Combined with our numerical results here, it may 
he concluded that the (xy,xy) scheme should not be taken 
in a general Gutzwiller calculation^ This conclusion is 
particularly important since the (xy,xy) scheme along 
Pier = is precisely the implicit scheme used by several 
authorsJ^^ii^ Among the remaining schemes, the (z,xy) 
and the diagrammatic schemes are probably preferred to 
the (z,z) scheme, since they provide better compromises 
between the accuracy in the small magnetization region 
and the large magnetization region. 

We have also performed a similar comparison between 
the Gutzwiller and the VMC results at a higher dopping 
of (5 = 0.125, where the optimal magnetization for the t-J 
model decreases to approximately 0.40. Most features of 
data discussed above in the 5 — 0.025 case continue to 
hold for the 5 ~ 0.125 case. The notable differences in the 
5 = 0.125 case are: 1. the erroneous minimum of {c\^Cj„) 
as function of (m) in the Gutzwiller approximation dis- 
appears for the pia = ni„ path; 2. the VMC results show 
a slightly stronger dependence on the choice of partic- 
ular paths in the parameter space; 3. Results from the 
Gutzwiller approximation along the paths = 1 and 
ijr = 1.1 significantly underestimated the optimal mag- 
netization for (H) (even though the shape of the curves 
still resemble that of the VMC results). The second and 
third differences are possibly caused by the flatness of 
the minimum, which is an indication that this dopping is 
close to the critical value of antifcrromagnet-paramagnet 
transition. 

It should be remarked that the VMC calculation along 
the ?/r = 1 line has been done previously by Yokoyama 
and Shiba^ and our results are consistent with theirs. It 
should also be noted that we have restricted our consid- 
eration to the antiferromagnet states in order to investi- 
gate the accuracy of the Gutzwiller approximation and its 
various reparation schemes in a clear and simple inhomo- 
geneous case. The antiferromagnctic state is somewhat 
artificial for the t-model at the dopping we considered, 
since the Nagaoka effeclj^^ suggests that the ground state 
should be ferromagnetic at low dopping, and indeed it 
can be checked that the ferromagnetic states have a lower 
expectation {c\„Cja) for both 6 = 0.025 and 5 = 0.125i^ 



IV. SUMMARY 

In this paper we have considered generalizing the 
Gutzwiller approximation to the case of an inhomogc- 
neous system and have found it useful to introduce extra 
spin-and-site-dependent fugacity factors. Wc derived the 
corresponding Gutzwiller factors from a configuration- 
counting approach in the appendix and explained its 
physical intuitions in the main text. The inclusion of 
fugacity factors reconcile the seemingly contradictory 
choices of Gutzwiller factors in the literature. Specifi- 
cally, different Gutzwiller factors that appear in the liter- 
ature correspond to different implicit choices of fugacity 



factors. This fact is particularly important when compar- 
ing results from the Gutzwiller approximation to those 
from other approaches, such as variational Monte-Carlo. 
The derivation and discussion of the Gutzwiller factors 
also show that the Gutzwiller approximation generally 
breaks the rotational symmetry of the trial wavcfunc- 
tion. Specifically, different components of the spin-spin 
interaction are rcnormalized differently and remain dif- 
ferent even in the homogeneous case, contrary to the as- 
sertions from the bulk of the literature. We proposed 
several possible schemes to remedy this defect of the 
Gutzwiller approximation, some of which are implicitly 
applied by various authors in the literature. It should be 
noted that these remedies are, strictly speaking, outside 
the scope of the configuration-counting approach of the 
Gutzwiller approximation. To compare the accuracy of 
the Gutzwiller approximation from various choices of fu- 
gacity factors and various choices of reparation schemes, 
we perform calculations for the two-dimensional square- 
lattice antiferromagnet and compare with the results 
from variational Monte-Carlo. Stated in general terms, 
the "lessons" learnt from the comparison are: 

1. There are regions in parameter space where the 
Gutzwiller approximation is erroneous. In particu- 
lar, one should avoid the parameter space where the 
pre-projected wavefunction and the fugacity factors 
each alone produces large inhomogeneity, but which 
the two counteract each other strongly. 

2. In general, one cannot assume gjxy and gjz (de- 
fined in Eq. [SHI31) to be equal, which otherwise 
can lead to qualitatively and quantitatively erro- 
neous results. It may be advisable, however, to 
take (7j| = gji = gjxy based on the consideration 
of rotational invariance. 

3. The energy expectation in both the t model and t- 
J model seems to depend strongly on physical pa- 
rameters, namely the single-site spin-specific den- 
sity {fiia), and weakly on whether such density is 
produced from the fugacity yia or from the pre- 
projected wavefunction I'i/'o)- States that produce 
the same {fiia) with different yia thus provide a 
possible consistency check. 

In short, it is advisable to use the (z,xy) scheme de- 
fined above Eq. [25] in calculations, and perform the cal- 
culations along at least two different paths in parameter 
space (conveniently the yr = 1 and {fii^) = {fiia)o paths) 
for a consistency check. 
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APPENDIX A: EXACT EXPRESSIONS FOR 
EXPECTATIONS OF GUTZWILLER 
PROJECTED WAVEFUNCTIONS 

We shall for simplicity assume the pre-projected wave- 
function |'0o) to be normalized, canonical and spin- 
definite. Explicitly, let |0) denotes the vacuum state, 
and let aka be linear combinations of Cja such that 



(0|afc4,|0) = Skk'- 
form: 



Then, \ipo) is assumed to take the 



i^o)-n4'Tn4j0)> 



(Al) 



where Egj^ E^j and E^j are projection operators onto 
empty j-site, up-spin j-site, and down-spin j-site respec- 
tively. Note that these projection operators are orthog- 
onal to each other (i.e., EajEpj = Sa/jEfjj), and that 



It is then an easy exercise to show that Gj = 
yio-"^lcr(^iff^]^), and hence for distinct sites ij, 



We shall also assume that there is some superlattice 
{R} such that {hj^^)Q = {nj+u,a)o and yj^a = yj+R.cr- 
Now, define Gj on site-j as: 

= (1 - njt)(l - fiji) + %t^jt(1 - %-J.) + - 



GjCj^Cja-Gj — UjaEuj 
GjGiA^CiaGiGj ~ yjayicr{Cjac]jg){cl„Cia-){cig-cl^ J 



GjGiC^j^Cj^clgCiaGiG- 



(A3) 



Eoj +2/jT-^b- +yjl^U 



(A2) 



From the orthogonality of the projection operators in 
Eq. lA2l and the explicit form of wavefunction assumed in 
Eq. lAll we have: 



im = (V'oi n^'^iv'o) = E (v^oi n yit^T. n yh^uU ^^M'o) 
= E (uyn] (iiy'] (v^(iin^^Tn(i-"^T)iV'^)(V',iin^^-in(i-"^i)i^o) 

{A,B) \jeA J \jGB / jeA ](A j£B j(B 

Here \iPq) = Ilfc '^Icrl'^)' with dj^^ defined as in Eq. lAU and the sum is over all subsets A,B of the set of lattice sites 
C such that 1^1 = E,(^o|",tI^o> = ^T> 1^1 = E,(^ol^jilV'o> = N^, and ^nS = 0. 
Taking also Eg. lASI into account, we have, analogously. 



(AM) \jeA J \jeB J jeA j^A jeB j(f:B 

i£A 



{A,B) \eeA\{t} J \eeB I i^A\{i] £^^u{j} fee i^B 

i^A.j<tAvjB 

(A6) 

= E ( n ) ( n ) (v-^i n n(i - ^n)i^(S>(^di n n(i - "n)iV'o^> (at) 

(AB) \^e-4 / XlaB / li^A l(^A IdB li_B 

i,j<^A 
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E 

{A.B) 




n (1 - "^t)i^(1)(^oi n 11(1 - n.^M^) 



(A8) 



(AM) 
i£Aj£B 




e&B\{j} e<^Bu{i} 



r 



(A9) 



and which the results for {iplc'j ^Cii\tp) , {ip\cl^Ciic]j^Cji\ip) 

and {iplcl-^CiiCj^Cjil^) can be obtained from the above 
with the substitution t*^! and A ^ B. 

It is known^ that the expectations appearing in 
Eg. IA4| - |A9l can be written in terms of determinants with 
elements of the form {ipolclcjltpo) ■ To summarize, 

(V'o I n n (1 - ^^-)l^o ) =det (AlO) 
ees e^s 

(V-o Isl^^^Il"^- 11(1 - ri,.)|0 = det4^) (All) 
ies\{i} etsu{j} 



where in Eq. lAllI it is assumed that i £ S while j ^ 
S, and that in both equations it is assumed that the 
constraint \S\ = X^iiV'ol^-io-IV'o) is met. Usa is an |£|-by- 
|£| matrix, whose entries are given by: 



when i € S 



-^jj - (V'oICcjtIV'o) whcni^5 



(A12) 



(ii) 

The matrix Ug^ is obtained from Usa by first exchang- 
ing the i-th and j-th column of Ug^ and then removing 
the j-th row and the j-th column from the resulting ma- 
trix. 



is expressed in terms of creation and annihilation opera- 
tors only on sites X = {ii, i„}, we first expand the ex- 
pectation in configuration basis as in Eq. IA4HA9I Then, 
each individual term in the expansion would depend on 
the details of the configuration. The Gutzwiller approxi- 
mation amounts to neglecting the non- combinatorial de- 
pendence on sites other than {«!,...,«„} in the quantum- 
mechanical part of each term. Hence it can be viewed as 
an approximation in evaluating the determinants det Usa 
and det C/l*i'' . For a canonical wavefunction defined on 

OCT 

a superlatticc, by considering the permutation expansion 
of the determinants, it can be seen that in general the 
off-diagonal terms depend on the correlations in the con- 
figuration while the purely diagonal term depends only 
combinatorially on the configuration. In other words, if 
S and S' contain the same number of sites on each sub- 
lattice, then the purely diagonal term in the permutation 
expansion of det Usa and det Us'a- agrees, while the off- 
diagonal terms in general disagree. The case for det U^^ 
is analogous. Hence, in the Gutzwiller approximation, we 
approximate: 



det Usa « det Usa-x x Y[ [Usa]j 



dct4^')«dct[/(^;jXn[C/54 

3ii 



(Bl) 



APPENDIX B: GUTZWILLER 
APPROXIMATION AND THERMODYNAMIC 
LIMIT 

The precise content of the Gutzwiller approximation 
can now be stated. When evaluating {ip\Q\ip), where Q 



where Uscr;X is the n-hy-n subblock of Usa whose ele- 
ments are those in Usa that connects one lattice sites in 
X to another, and that U^^^i is obtained from Usa;i in 
the same way which det U^ J is obtained from Usa- 
Applying the approximation. Eg . I A4I - I A9I simplifies to: 



J 



E 

(A,B)3<^A 



n n (1 " '^Ji) n yn^^oi ^^oi) 

liA jeB j(^B 



(B2) 



(V'ici|C,ti?/') E n 2^jt"jt n (1 " "jt) n yn^^^i n(i ~ "^'^^ 

3<fA 3&B j<^B 



(B3) 



{A,B}j€A 
ieA 
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(V'lcjt^iTl 



iA,B) 



£eA\{i} i^AuU} 



(B4) 



{A,B) 
i,j£A 



eGA\{u3} 



eeB 



l^B 



(B5) 



(-4,6) 



fe^\{i} 



(B6) 



(B7) 



where 



in above and henceforth. 



With our assumption of a supcrlattice structure, the lattice can be divided into k sublattices. Let M = I^CI/k be 
the number of unit supcrcell, then a httlc bit of combinatorics (c.f. Eg . [T51 and [TB|) gives: 



(^1^) = E n 



M! 



(B8) 



= E 



F{{ai„}) (B13) 



where / label the sublattices, and F{{ai„}) is defined 
in the obvious way. The sum is over all non-negative 
sets of {a/|} and {a/j,} such that X^/'^^T ~ -^T ^^"^ 
Y^iO-ii = N^. In a similar manner, it can be checked The corresponding expressions when P and Q are iden- 
that Eq. [BSHBTI can be written in the general form tical are similar, and yield the sa me formu la once the 



Y^l^j^yi- ■ .)F{{aicr}). Specifically, for site-i is on sub- thermodynamic hmit is taken (Eq. |B14|[BT3 below). 



lattice P and site-j is on sublattice Q, with P and Q 
distinct 



{aic} 



As discussed in the main text, in the thermodynamic 
limit, whereby |£| —>■ oo, the sum ^^{a^ } '^^^ tie replaced 
by the single term in which F{{aj^}) attains maximum 
(B9) under the constraints "^jCli^ ~ No-. The maximization 
condition results in Eq. [18] in the main text, with p/o- 
formally defined as pj„ = aja-/M in the present context. 
Substituting the pia into Eq. IBSffBlSI then gives: 



ap„{M - OQi - aQi)ypayQa 



{flier) = PPc 



(B14) 



E M^{yl^np,)il^nQ,) 



Fi{ai,}) 



(BIO) 



^ a'PlaQ^yp^yplyQ^yQlFi{aI<T}) 

W }^^^d - np^)iyl^npi){l - nQ|)(y^^nQ|) 

(Bll) 



(^L^«0"^]<T''^iO'') 



yp„ npa{l-nQ^) 



(B15) 
(B16) 



PPa PQa' /4 „ 4 . , 

npcrnQa' 

yPayQcr npa{l ~ np<j)nQcr(l - riQff) 

(B17) 



'\c\„Ciac]„Cja\^) = (cLci<rc]^Cj<,)o E 



{aia} 



«p ^({a/ }) ^^^'^^ site-i is on sublattice P while site-j is on sub- 
lattice Q (P and Q may be identical). Using Eq. [T51 
it can be checked explicitly that the hermiticity of ex- 



JVPnp^UQa 



(B12) pectation values is preserved under the approximation. 
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Specifically, the R.H.S. of Eq. IB 151 and Eq. IB 171 remain the the yi„ in these equations. Upon simplifications and 
unchanged upon P ^ Q. This allows us to eliminate rearrangements this yields Eg . in the main text. 
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